Last week I was on Gran Canaria for a vacation We took hiking trips and explored the island by car and on foot. So, what better way to keep up the holiday spirit a while longer than to visualize all the places we went!?
I am combining the location data collected by
I am using ggplot and ggmap. I centered the map to the midpoint of Gran Canaria as given by http://www.travelmath.com/island/Gran+Canaria, which is 27° 57’ 54" N / 15° 35’ 59" W. Converting to decimal with http://www.rapidtables.com/convert/number/degrees-minutes-seconds-to-degrees.htm that is 27.965 & 15.59972.
library(ggplot2)
library(ggmap)
map_theme <- list(theme(legend.position = "top",
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_rect(fill = "white"),
panel.border = element_blank(),
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size = 18)))
map <- get_map(c(lon = -15.59972, lat = 27.965), zoom = 10, maptype = "terrain-background", source = "google")
#map <- get_map(c(lon = -15.59972, lat = 27.965), zoom = 10)
The GPS tracks in .gpx format were downloaded from the device (Garmin) the same way as I had already done with the tracks from our US/Canada roadtrip last year. Garmin’s .gpx tracks were not in a standard format that could be read with readGPX(), so I had to use htmlTreeParse() from the XML library.
library(XML)
myfiles <- list.files(path = "Takeout_2017_March", full.names = TRUE)
myfiles <- myfiles[grep(".gpx", myfiles)]
for (i in 1:length(myfiles)) {
tryCatch({
pfile <- htmlTreeParse(myfiles[i], useInternalNodes = T)
}, error = function(e){cat("ERROR\n")})
if (exists("pfile")) {
# Get all elevations, times and coordinates via the respective xpath
elevations <- as.numeric(as.character(xpathSApply(pfile, path = "//trkpt/ele", xmlValue)))
times <- xpathSApply(pfile, path = "//trkpt/time", xmlValue)
coords <- xpathSApply(pfile, path = "//trkpt", xmlAttrs)
# Extract latitude and longitude from the coordinates
lats <- as.numeric(as.character(coords["lat",]))
lons <- as.numeric(as.character(coords["lon",]))
# Put everything in a dataframe and get rid of old variables
geodf <- data.frame(lat = lats, lon = lons, ele = elevations, time = times)
if (i == 1) {
geodata <- geodf
} else {
geodata <- rbind(geodata, geodf)
}
}
rm(pfile)
}
geodata_all <- geodata
# Transforming the time column
geodata_all$time <- as.character(strptime(geodata_all$time, format = "%Y-%m-%dT%H:%M:%SZ"))
# ordering by date
library(dplyr)
geodata_all <- arrange(geodata_all, time)
geodata_gc <- filter(geodata_all, lat < 28.5 & lat > 27.5) %>%
filter(lon > -16 & lon < -15)
The hiking tracks we followed came mostly from the Rother Wanderführer, 7th edition from 2016. They were in standard .gpx format and could be read with readGPX().
library(plotKML)
myfiles <- list.files(path = "hiking", full.names = TRUE)
myfiles <- myfiles[grep(".gpx", myfiles)]
for (i in 1:length(myfiles)) {
t <- readGPX(myfiles[i])
geodf <- data.frame(lon = t$tracks[[1]][[1]]$lon,
lat = t$tracks[[1]][[1]]$lat,
ele = t$tracks[[1]][[1]]$ele,
tour = names(t$tracks[[1]]))
if (i == 1) {
geodata <- geodf
} else {
geodata <- rbind(geodata, geodf)
}
rm(pfile)
}
geodata_tracks <- geodata
The Google location data I loaded the same way as in my post about plotting your Google location history.
library(jsonlite)
system.time(x <- fromJSON("Standortverlauf.json"))
# extracting the locations dataframe
loc = x$locations
# converting time column from posix milliseconds into a readable time scale
loc$time = as.POSIXct(as.numeric(x$locations$timestampMs)/1000, origin = "1970-01-01")
# converting longitude and latitude from E7 to GPS coordinates
loc$lat = loc$latitudeE7 / 1e7
loc$lon = loc$longitudeE7 / 1e7
# keep only data from Gran Canaria
loc_gc <- filter(loc, lat < 28.5 & lat > 27.5) %>%
filter(lon > -16 & lon < -15)
All three datasets had information about the latitude/longitude coordinate pairs and of the elevation of each observation, so I can combine these three attributes from the three location data sources.
geodata_tracks$ele <- as.numeric(as.character(geodata_tracks$ele))
data_combined <- data.frame(lat = c(geodata_gc$lat, loc_gc$lat, geodata_tracks$lat),
lon = c(geodata_gc$lon, loc_gc$lon, geodata_tracks$lon),
ele = c(geodata_gc$ele, loc_gc$altitude, geodata_tracks$ele),
track = c(rep("GPS", nrow(geodata_gc)), rep("Google", nrow(loc_gc)), rep("Hiking", nrow(geodata_tracks))))
data_combined <- data_combined[!duplicated(data_combined), ]
ggmap(map) +
geom_point(data = data_combined, aes(x = lon, y = lat, color = track), alpha = 0.3) +
scale_color_brewer(palette = "Set1") +
map_theme
ggmap(map) +
geom_point(data = na.omit(data_combined), aes(x = lon, y = lat, color = ele)) +
scale_color_gradient2(low = "lightblue", mid = "blue", high = "red", name = "Elevation") +
map_theme
library(dplyr)
geodata_gc$time <- as.POSIXct(geodata_gc$time, format = "%Y-%m-%d %H:%M:%S", tz = "GMT")
data_combined_2 <- geodata_gc
# removing duplicate entries
data_combined_2 <- data_combined_2[!duplicated(data_combined_2), ]
# order by time
data_combined_2 <- data_combined_2[order(data_combined_2$time, decreasing = FALSE), ]
# Shifting vectors for latitude and longitude to include end position
shift.vec <- function(vec, shift) {
if (length(vec) <= abs(shift)) {
rep(NA ,length(vec))
} else {
if (shift >= 0) {
c(rep(NA, shift), vec[1:(length(vec) - shift)]) }
else {
c(vec[(abs(shift) + 1):length(vec)], rep(NA, abs(shift)))
}
}
}
data_combined_2$lat.p1 <- shift.vec(data_combined_2$lat, -1)
data_combined_2$lon.p1 <- shift.vec(data_combined_2$lon, -1)
# Calculating distances between points (in metres) with the function pointDistance from the 'raster' package.
library(raster)
data_combined_2$dist.to.prev <- apply(data_combined_2, 1, FUN = function(row) {
pointDistance(c(as.numeric(as.character(row["lat.p1"])),
as.numeric(as.character(row["lon.p1"]))),
c(as.numeric(as.character(row["lat"])), as.numeric(as.character(row["lon"]))),
lonlat = T) # Parameter 'lonlat' has to be TRUE!
})
round(sum(as.numeric(as.character(data_combined_2$dist.to.prev)), na.rm = TRUE) * 0.001, digits = 2)
## [1] 1009.6
library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
library(raster)
srtm <- getData("SRTM", lon=-15.59972, lat=27.965)
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
raster::persp(srtm,
xlab = "Longitude", ylab = "Latitude", zlab = "Elevation",
phi = 45, theta = 0,
expand = 0.1,
col = "forestgreen", border = NULL,
box = TRUE)
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
e <- extent(min(data_combined$lon) - 0.2, # xmin
max(data_combined$lon) + 0.1, # xmax
min(data_combined$lat) - 0.1, # ymin
max(data_combined$lat) + 0.1) # ymax
srtm_c <- crop(srtm, e)
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
raster::persp(srtm_c,
xlab = "Longitude", ylab = "Latitude", zlab = "Elevation",
phi = 45, theta = 45,
expand = 0.5,
col = "forestgreen", border = NULL,
box = TRUE)
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
raster::persp(srtm_c,
xlab = "Longitude", ylab = "Latitude", zlab = "Elevation",
phi = 45, theta = 225,
expand = 0.5,
col = "forestgreen", border = NULL,
box = TRUE)
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
x <- terrain(srtm_c, opt = c("slope", "aspect"), unit = "degrees")
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
plot(x)
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
slope <- terrain(srtm_c, opt = "slope")
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
aspect <- terrain(srtm_c, opt = "aspect")
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
hill <- hillShade(slope, aspect, 40, 270)
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
plot(hill, col = grey(0:100/100))
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
library(rgdal)
library(rasterVis)
library(rgl)
library(htmlwidgets)
options(rgl.printRglwidget = TRUE)
open3d()
## glX
## 1
plot3D(srtm)
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files